In class exericse 07

Time series

Frostbear https://sg.linkedin.com/in/farahfoo (SMU Masters in IT business (Fintech and Analytics))https://scis.smu.edu.sg/master-it-business
2022-02-26

TABLEAU EXERCISE

Map here

R EXERCISE

Ex 1:Mapping Geospatial Point Data with R

To load packages required today:

packages = c('tidyverse', 'sf', 'tmap')

for(p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p, character.only = T)
}

sgpools <- read_csv ("data/aspatial/SGPools_svy21.csv")

head (sgpools, 3)
# A tibble: 3 x 7
  NAME   ADDRESS POSTCODE XCOORD YCOORD `OUTLET TYPE` `Gp1Gp2 Winning~
  <chr>  <chr>      <dbl>  <dbl>  <dbl> <chr>                    <dbl>
1 Livew~ 2 Bayf~    18972 30842. 29599. Branch                       5
2 Livew~ 26 Sen~    98138 26704. 26526. Branch                      11
3 Sport~ Lotus ~   738078 20118. 44888. Branch                       0
sgpools_sf <- st_as_sf(
  sgpools,
  coords = c("XCOORD", "YCOORD"),
          crs=3414)

head (sgpools_sf, 3)
Simple feature collection with 3 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 20117.93 ymin: 26525.7 xmax: 30841.56 ymax: 44888.06
Projected CRS: SVY21 / Singapore TM
# A tibble: 3 x 6
  NAME       ADDRESS           POSTCODE `OUTLET TYPE` `Gp1Gp2 Winning~
  <chr>      <chr>                <dbl> <chr>                    <dbl>
1 Livewire ~ 2 Bayfront Avenu~    18972 Branch                       5
2 Livewire ~ 26 Sentosa Gatew~    98138 Branch                      11
3 SportsBuz~ Lotus Lounge, Le~   738078 Branch                       0
# ... with 1 more variable: geometry <POINT [m]>
tmap_mode("view")

tm_shape(sgpools_sf) +
tm_bubbles(col = "green",
            size = 0.5,
            border.col = "black",
            border.lwd = 1,
           interactive = TRUE)

Faceting 2 maps for synchronisation view

sync = True will ensure the 2 maps move in tandem

tm_shape(sgpools_sf) +
tm_bubbles(col = "OUTLET TYPE",
            size = "Gp1Gp2 Winnings",
            border.col = "black",
            border.lwd = 0.5,
           interactive = TRUE) +
  tm_facets(by= "OUTLET TYPE",
            nrow = 1,
            sync = TRUE)

Ex 2: Choropleth Mapping with R

use st_read function to convert the .shp file into 1 single tibble object, to get an administrative boundary map, and a column of geometry (string of coordinates)

mpsz <- st_read(dsn = "data/geospatial",
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\farahfoo\VA\In_Class_ex\In-class_Ex07\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
head (mpsz, 1)
Simple feature collection with 1 feature and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 30794.28 ymin: 28369.47 xmax: 32362.39 ymax: 30140.01
Projected CRS: SVY21
  OBJECTID SUBZONE_NO    SUBZONE_N SUBZONE_C CA_IND   PLN_AREA_N
1        1          1 MARINA SOUTH    MSSZ01      Y MARINA SOUTH
  PLN_AREA_C       REGION_N REGION_C          INC_CRC FMEL_UPD_D
1         MS CENTRAL REGION       CR 5ED7EB253F99252E 2014-12-05
    X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area
1 31595.84 29220.19   5267.381    1630379
                        geometry
1 MULTIPOLYGON (((31495.56 30...

Read in the population data

popagsex <- read_csv("data/aspatial/respopagsex2000to2018.csv")

head (popagsex, 2)
# A tibble: 2 x 6
  PA         SZ        AG     Sex       Pop  Time
  <chr>      <chr>     <chr>  <chr>   <dbl> <dbl>
1 Ang Mo Kio Cheng San 0_to_4 Males     810  2000
2 Ang Mo Kio Cheng San 0_to_4 Females   700  2000

prepare the population data for use

popagsex2018_male <- popagsex %>%
filter(Sex == "Males") %>%
filter(Time == 2018) %>%
spread(AG, Pop) %>%
mutate(YOUNG = `0_to_4`+`5_to_9`+`10_to_14`+
`15_to_19`+`20_to_24`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13]) + 
                          rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)/`ECONOMY ACTIVE`) %>%

# this converts the PA and SZ to upper case to left_join later on
mutate_at(.vars = vars(PA, SZ),.funs = funs(toupper)) %>%
  
select(`PA`, `SZ`, `YOUNG`,`ECONOMY ACTIVE`, `AGED`,`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` >0)

head (popagsex2018_male, 2)
# A tibble: 2 x 7
  PA         SZ          YOUNG `ECONOMY ACTIVE`  AGED TOTAL DEPENDENCY
  <chr>      <chr>       <dbl>            <dbl> <dbl> <dbl>      <dbl>
1 ANG MO KIO ANG MO KIO~   660             1290   310  2260      0.752
2 ANG MO KIO CHENG SAN    3110             8020  2440 13570      0.692

left join the 2 tables using SUBZONE_N and SZ

glimpse (popagsex2018_male)
Rows: 233
Columns: 7
$ PA               <chr> "ANG MO KIO", "ANG MO KIO", "ANG MO KIO", "~
$ SZ               <chr> "ANG MO KIO TOWN CENTRE", "CHENG SAN", "CHO~
$ YOUNG            <dbl> 660, 3110, 2900, 2370, 890, 1860, 1050, 235~
$ `ECONOMY ACTIVE` <dbl> 1290, 8020, 7610, 6370, 1730, 4800, 2200, 5~
$ AGED             <dbl> 310, 2440, 2710, 2130, 540, 1430, 640, 2020~
$ TOTAL            <dbl> 2260, 13570, 13220, 10870, 3160, 8090, 3890~
$ DEPENDENCY       <dbl> 0.7519380, 0.6920200, 0.7371879, 0.7064364,~
mpsz_agemale2018 <- left_join(mpsz,popagsex2018_male,by = c("SUBZONE_N"="SZ"))

head (mpsz_agemale2018, 2)
Simple feature collection with 2 features and 21 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 28160.23 ymin: 28369.47 xmax: 32362.39 ymax: 30177.73
Projected CRS: SVY21
  OBJECTID SUBZONE_NO    SUBZONE_N SUBZONE_C CA_IND   PLN_AREA_N
1        1          1 MARINA SOUTH    MSSZ01      Y MARINA SOUTH
2        2          1 PEARL'S HILL    OTSZ01      Y       OUTRAM
  PLN_AREA_C       REGION_N REGION_C          INC_CRC FMEL_UPD_D
1         MS CENTRAL REGION       CR 5ED7EB253F99252E 2014-12-05
2         OT CENTRAL REGION       CR 8C7149B9EB32EEFC 2014-12-05
    X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area     PA YOUNG ECONOMY ACTIVE
1 31595.84 29220.19   5267.381  1630379.3   <NA>    NA             NA
2 28679.06 29782.05   3506.107   559816.2 OUTRAM   620           2210
  AGED TOTAL DEPENDENCY                       geometry
1   NA    NA         NA MULTIPOLYGON (((31495.56 30...
2 1160  3990  0.8054299 MULTIPOLYGON (((29092.28 30...

note that this map has missing values as there are some ares with no population

tmap_mode("plot")

qtm(mpsz_agemale2018,fill ="DEPENDENCY")

Drawing the map with tm_fill and tm_border

tm_shape(mpsz_agemale2018) +
tm_polygons()
tm_shape(mpsz_agemale2018) +
tm_polygons("DEPENDENCY")

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY")
tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY")+
  tm_borders(lwd = 0.1,
             alpha = 1)

Data classification methods of tmap

Most choropleth maps employ some method of dataclassification. The point of classification is to take alarge number of observations and group them into data ranges or classes.

Compare the 2 styles below: quantile vs equal, the intervals are different

tm_shape(mpsz_agemale2018) + 
  tm_fill("DEPENDENCY",n =8,style ="quantile") +
tm_borders(alpha =0.5) +
  tm_layout(main.title ="Quantile interval distribution")

tm_shape(mpsz_agemale2018) + 
  tm_fill("DEPENDENCY",n =8,style ="equal") +
tm_borders(alpha =0.5) +
  tm_layout(main.title ="Equal interval distribution")

adding colour palette

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY", n = 6,style ="quantile",palette ="Blues") +
  tm_borders(alpha =0.5) 

just add a “-” sign infront to inverse the colour scale
tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY", n = 6,style ="quantile",palette ="-Blues") +
  tm_borders(alpha =0.5) +
  tm_layout(main.title ="Reversing colour scheme",
  main.title.position ="center",
main.title.size =1,
legend.height =0.45,
legend.width =0.35,
legend.outside =FALSE,
legend.position = c("right", "bottom"),
frame =FALSE)

Add ing Cartographic Furniture into the tmap

Add in compass, scalebar, gridlines into the map

tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",style ="quantile",palette ="Blues",title ="No. of persons") +
tm_layout(main.title ="Distribution of Dependency Ratio \nby planning subzone",
main.title.position ="center",main.title.size =1.2,
legend.height =0.45,
legend.width =0.35,
frame =TRUE) +
tm_borders(alpha =0.5) +
tm_compass(type="8star", size =2) +
tm_scale_bar(width =0.15) +
tm_grid(lwd =0.1, alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from Urban Redevelopment Authorithy (URA)\n and Population data from Department of Statistics DOS",
             position = c("left","bottom"))

Faceting the maps and custom view to show applicable area only

tm_shape(mpsz_agemale2018) +
  tm_fill("DEPENDENCY",
          style = "quantile",
          palette = "Blues",
          thres.poly = 0) + 
  tm_facets(by="REGION_N",
            free.coords=TRUE,
            drop.shapes=TRUE) +
  tm_layout(legend.show = FALSE,
            title.position = c("center",
                               "center"), 
            title.size = 20) +
  tm_borders(alpha = 0.5)